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LARGE EDDY INTERACTIONS IN A TURBULENT CHANNEL FLOW 


S. K. Hong* 

Ames Research Center 


SUMMARY 


The dynamic processes of large eddies in a turbulent channel flow have been 
examined by utilizing an orthogonal expansion of the velocity fluctuation, known in 
the literature as the Proper Orthogonal Decomposition Theorem. The mathematical 
form of these functions is unknown, in contrast to the Fourier analysis. Attention 
is focused on the nonlinear, turbulence-turbulence interaction process in the dynam- 
ical equation for large eddies (the first term in the expansion). The nonlinear 
interactions of the components of the first mode are treated exactly, but influences 
of higher modes are modeled. This requires adjustment of both the skewness and the 
effective Reynolds number so that the energy equilibrium of the large eddies is 
ensured when the mean velocity distribution is assumed known from experiments. 
Computational results show that the first mode contributes significantly to turbu- 
lent intensities and possesses a structural and statistical character similar to 
that of the entire flow. 


INTRODUCTION 


Important in the engineering predictions of inhomogeneous, turbulent shear 
flows are the mean velocities and the local structures of turbulence under given 
initial and boundary conditions. The most popular prediction methods first assume 
that the Navier-Stokes equations are adequate for describing turbulent flow on an 
instantaneous basis and then proceed to develop statistical equations for the vari- 
ous turbulent moments, including the Reynolds stresses. However, these equations 
involve more turbulent moments than equations that exist for them, forming an open 
system. Various levels of closure schemes have been proposed: the zero-equation 

model (Cebeci and Smith, 1974), the two-equation model (Saffman, 1970; Jones and 
Launder, 1972), and the Reynolds stress equation model (Hanjalic and Launder, 1972; 
Mellor and Herring, 1973). Each of the foregoing methods requires the introduction 
of several empirical constants with respect to various turbulent processes and 
provides only approximate predictions of the nature of individual turbulent pro- 
cesses arising in a given flow. An approach that examines the dynamics of the 
turbulence may require less reliance on modeling. Earlier, the author used the 
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Proper Orthogonal Decomposition Theorem (PODT) (Loeve, I960; Lumley, 1967) to 
develop a framework for predicting turbulence structural quantities by modeling 
nonlinear eddy-eddy interaction terms. At the same time, the framework was used to 
study the implications of the modeling adopted. The framework, called the Large 
Eddy Interaction Model (LEIM) (Hong, 1983), is based on the concept that a few 
properly identified modes in the PODT may be used to represent the technologically 
important quantities of the turbulent flow field. 

It should be noted that Lumley's (1967) primary purpose was to define, unambig- 
uously, the meaning of an "eddy." Given two-point velocity correlations from exper- 
iments, the PODT has been applied to pipe flow (Bakewell and Lumley, 1967), a wake 
(Payne and Lumley, 1967), and a flat plate boundary layer (Lemmerman and Payne, 

1977) as a means of extracting the features of a dominant eddy. In a similar 
approach using the results from a computationally simulated fully developed channel 
flow, Moin (1984) specifically investigated the number of modes necessary in the 
PODT to reproduce the turbulent intensities and shear stress. Moin (1984) shows 
that, m the case of shear stress, it takes the first 15 terms in the PODT represen- 
tation before the calculated stress distribution matches the value simulated earlier 
(Moin and Kim, 1982) across the boundary layer. However, when only the wall region 
is examined, the sum of the first three modes yields the "experimental" results 
quite well. This small number of modes needed to describe the turbulence gave the 
author further encouragement to extend the PODT approach, as was done in the LEIM, 
to be a predictive tool. 

The Large Eddy Interaction Model has been applied m the past to variously 
curved wall boundary layer flows (Hong and Murthy, 1983, 1984a, 1984b). The frame- 
work has proved to be useful, among other features, in establishing directly the 
manner in which (a) anisotropy can arise and change and (b) turbulent transport is 
affected by the addition or removal of an extra strain m those complex flows. The 
basic procedure of the LEIM consists of the following steps: (a) decomposing the 

velocity fluctuations into orthogonal functions with random coefficients, (b) con- 
structing dynamical equations for those functions, (c) identifying the first mode as 
an organized structure that contributes most to the energy (Lumley, 1967, 1981), and 
(d) evaluating the large eddy which interacts with the mean flow and the eddy-eddy 
interactions. However, all the nonlinear terms m the LEIM were modeled in a linear 
form utilizing either an anisotropic eddy viscosity or a diffusion velocity. In 
this process, three empirical constants were introduced in the closure and were then 
determined by matching shapes between normalized Reynolds stresses, calculated from 
the first mode, and experimental measurements. In view of the emphasis of past 
applications of the LEIM on evaluating the normalized structure of the first mode, 
it was primarily a diagnostic method. 

In the present work to develop a predictive method, the turbulent transport 
processes have been re-examined and retained in their nonlinear form. This mini- 
mizes the dependence on turbulence modeling and allows evaluation of the magnitudes 
of the moments. The applicability of the new transport model has been illustrated 
in a channel flow that is inhomogeneous in the direction normal to the wall. How- 
ever, the computed results shown here are restricted to the use of only a single 
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mode in the decomposition. As stated earlier, this forced the introduction of the 
skewness and the effective Reynolds number as parameters of the problem. Although 
the current work has not yet demonstrated the uniqueness of a set of these param- 
eters, results based on different sets of these parameters (satisfying the energy 
equilibrium of large eddies) yield Reynolds stresses that are less than a few per- 
cent apart over the entire channel. Structural and statistical properties of the 
large eddies are then studied from the solution of the first mode. Findings are 
consistent with the earlier studies on the use of the PODT (that the first mode 
exhibits the structural nature of the averaged turbulence moments and the statisti- 
cal nature of the random turbulence). The long term objective is that the LEIM 
framework be developed further, not only to gain the detailed phenomenological 
insight into turbulence, but also to provide a means to predict turbulent quantities 
of engineering interest. 

The author appreciates the suggestions provided by Dr. M. W. Rubesin during the 
course of the work, and his comments on the manuscript. 


LARGE EDDIES AND THEIR INTERACTIONS 


The mathematical definition of a turbulent eddy as proposed by Lumley (1967), 
and the method of relating the structure of turbulence to that of a large eddy as 
developed by Lumley (1967, 1981) and Townsend (1976, 1980), are applied here to 
turbulent shear flow that is inhomogeneous at least in one spatial direction. A 
closure assumption is introduced for the nonlinear eddy-eddy interaction term aris- 
ing in the dynamical equation of the large eddies. 


Proper Orthogonal Decomposition Theorem 

One can consider a decomposition of the instantaneous velocity, 0 ^ , into a mean 
value and a fluctuation as 



+ u. 
i 


(D 


Here we assume that is an ensemble average of 0.. Then it is possible to 

write the velocity fluctuation in terms of orthogonal functions, 


{ *i 


n = 1,2,3,...}. That is, 


u^x.t) = a n <j/ n ^(x,t) 

n=1 


dx dt = 0 


(p * q; l = 1, 2, or 3) (2) 
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where a^, a^, • •• are random coefficients with units of velocity and uncorre- 

lated from one another, which means 


»„ * 0 

Vi = x< " )5 mnj 

6 is the Kronecker delta function defined as 
mn 


mn 


f 1 

to i 


if m = n 
if m * 0 


(3) 


The overbar, (~), represents the expectation of the event, ( ) , in the probability 
theory and the ensemble average in the turbulence theory. Equation (2) is a gener- 
alized Fourier series, accounting for the inhomogeneity of flows by means of dis- 
crete functions. It is assumed by definition that a^, a^, a ... are ordered 
such that 


x (1 > > x< 2) > x (3) 


> 0 


If one assumes U L in eq. (1) to be a time-averaged velocity, one needs a slightly 
different expansion from eq. (2) as described m appendix A. Equation (3) implies 
that the coefficients are completely random so that the correlation between the same 
coefficients is perfect, and between different coefficients is zero. In addition 
the <}>p's are assumed to be orthonormal functions, that is 


/‘i P M q> d * dt ■ 5 pq 


(4) 


It is noted that the orthonormality condition is imposed only for the u-component 
functions, and the limits of the integral in eq. (4) are quite arbitrary. The 
orthogonality condition also implies that none of <j>' n ^'s is identically zero 
(Tolstov, 1962, p. 41). To calculate the coefficients, a^, both sides of eq. (2) 
when 1 = 1 are multiplied by <j>!j n , and the resulting equation is integrated 
taking account of the orthonormality condition, eq. (4). The Fourier coefficients 
are then given by the following form. 


a 

n 


J* u(x,t)<t>!j n ) 


(x,t)dx dt 


(5) 



It can be shown, utilizing the randomness in a n and the orthogonality m 
that 4 / is related to the two-pomt velocity correlation, R^j, as 
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( 6 ) 


R i j(x l x' ;t,t' ) = u i (x,t)Uj(x' ,t' ) 


. i 

n=1 


(x,t)<t>j n ^ 




and 


J'r.^XjX' ;t,t' )4>^(x' ,t' )dx * dt' = x^ n ^ n ^(x,t) (7) 

where 4>^ n ^ are eigenfunctions with X^ as their eigenvalues of the system. 
Further, if one considers the integral of the turbulent kinetic energy from the 
u-component over the whole flow field, then 


/ i u 2 d5 dt = / 1 £ 


(n). 


(n) 

1 




dx dt 




( 8 ) 


(n) 


represents the kinetic energy content of the entire flow associated 
-mode. The nth mode distributes the energy, X^ , in space and time 


Thus, 
with 

according to its functional form. 


Formulation of Large Eddy Interaction Model 


In incompressible turbulent flow, decomposition of the Navier-Stokes equation 
into the mean and fluctuating parts leads to the following dynamical equation for 
the velocity fluctuation. 


3u . 3u. 3U 

1 .. l l 

3t + j 3x j + 3Xj U j 



(u i u j 


- Vj> 


1 3P 

P 3x i 



(9) 


where p is the pressure fluctuation. Introducing the orthogonal decomposition of 
the velocity fluctuation, eq. (2), into eq. (9), a dynamical equation is obtained 
for the eigenmode, , by the following procedure. First, u 1 is replaced by its 
expansion; secondly, both sides of the resulting equation are multiplied by a random 
coefficient, c* m ; thirdly, an average is obtained by utilizing the re latio n , 

a m a n = X^ n ^6 m ; and finally, the remaining equation is divided by ^/x ^ n ^ , which 
leads to the m ?ollowing as first given by Lumley (1967). 
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a{! n) ai (n) 3U, , , 

_i . „ 1 . _i r(n) 

3t °j 3Xj 3Xj * j 
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J |.p=i q=i 


a a a 
n P q 


3 X j 1 p =1 q =1 (X (n) X ( P ) X (< ’ ) ) 1/2 


_2~(n) 

3ir " 3 *i 

+ v ~ — 
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3X 


j 




(n = 1, 2, 3, ...) 


(10) 


where 


and 



(ID 


( n ) 

it ' is also the nth mode in the series expansion of the pressure fluctuation, 
that is, 


1 (1) (2) (3) 

p i 2 3 

It is noted that no constraints are imposed on other than determinacy. 

Equation (10) thus governs both the shape and intensity of each 4^-mode. The 
parameter, n, represents the order of modes where it is assumed that the first mode, 
n = 1, accommodates most of the energy in and of structure in the dimensionless 

function, <t> ^ . Similarly, the second mode, n = 2, is supposed to take up most of 
what is left, that is, u - ct^ . The index, l, may have 1, 2, or 3 corresponding 
to the streamwise direction (x), the local normal to the wall (y), or the spanwise 
direction (z), respectively. While it is of great interest to solve such a system 
of equations in general, attention is focused here on the lowest mode for which the 
dynamical equation becomes the following. 


3<J) 


( 1 ) 


3 4> 


( 1 ) 


3t 


+ U 


j 3X, 


3U. 
1 

3X 


r(1> 

*3 


arft 

j [p-- 


a.aa 

i p q 


q =1 <x ( 'V p) x (q V /2 


;^ q) } 


3 ; (1 > 

3Xi 


2r( 1 ) 


3 <t> 


+ V 


( 12 ) 


3x 
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The eddy-eddy interaction term in the above may be said to consist of interactions 
of large eddy with large eddy, of large eddy with small eddy, and of small eddies 
with themselves. 
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For the first eigenmode of pressure fluctuation appearing in eq. (12), one can 
develop the following Poisson's equation. 


a 2 ; (1) 

3X^ 


3U. 3<t> 


r(1) 


3X, 


3 X 


2 ^ CD CD 

3x73x7 ] 23 23 

k j l_p=1 q = 1 


a . act 

i p q 


r(q)r(p) 

h ( 1 n) 1 (p) l (q) ) , /a »J *k 




(13) 


On the right-hand side of the Poisson's equation are contributions due to interac- 
tion between mean velocity gradient and large eddy velocity gradient, and due to 
interactions among eddies of all sizes. Or, alternatively to eq. (13), continuity 
relation can be used in the form 


axj 


= 0 


(14) 


Equation (12) with eq. (13) or (14) forms a system of equations for the first 
mode, <|> . The nonlinear eddy interaction term in those equations involves higher 


1 ( 2 ) 


'I 
♦i 
or 


(3) 
>i 


(4) 
; i 


Thus in order to close the system of equations for 


m^es, <)k , , 

, it is necessary either to drop the interactions associated with higher modes 
to model the contributions from them. 


Transport Process 


The nonlinear interaction between the first mode is retained in its original 
form but the interactions involving modes higher than <{7 are modeled. A simple 
way of accounting for the effect of the higher modes is to group them together and 
to relate this effect to a known quantity. An eddy viscosity is introduced for this 
purpose (Townsend, 1976, 1980). 


Y* Y 1_B_3 


1 rOMD 

(1)3/2 *i *J 



(15) 

which has also been suggested by Lumley (1967). In eq. (15), v^, denotes an eddy 
viscosity. In the present analysis, is kept equal to a constant which is inde- 
pendent of mesh dimensions or the distance from the wall. This has an effect of 
reducing the effective Reynolds number by a factor of v/(v + v^). Although 
improved models (with variation m the y-direction) may be required, for simplicity 
we will use constant to achieve a steady state solution. 

U pon su bstituting eq. (15) into (12), one obtains a closed system of equations 
for j/^>. 
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( 16 ) 


3_ 

at 


id) 

*1 


_a 

ax 
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; ( 1 } 
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au. 

l 

3X , 




a 2 r( 1 ) 

3 ~(D , > 3 *i 

= _ ^ + ( V + v T ) 3~ i7T 

l j J 

“ 3 ~2 1/2 

where S (= 0 ^/( 0 ^-^ ) is the skewness factor of the random coefficient, c^. 

Assuming that the mean velocity is given, the system of equations involves then a 
structural parameter, S, and a stability parameter, v T , which need to be chosen. In 
order for the LEIM to be a completely predictive scheme, it requires incorporation 
of the mean momentum equation for along with eqs. (14) and (16). For the 

present, however, emphasis will be placed on how the large eddies interact and react 
to a known mean flow field. 


APPLICATIONS TO CHANNEL FLOW 


Fully developed turbulent channel flow has acquired a large data base over the 
years. Accordingly, two-dimensional channel flow has been chosen to demonstrate the 
validity of the approximation for nonlinear, eddy-eddy interaction terms, as pro- 
posed in eq. (15), and to establish the contribution of the first mode to various 
statistical turbulence quantities. In the fully developed region of the channel, 
the mean velocity is one-dimensional. It is dependent only on the normal coordi- 
nate, y, where y = 0 corresponds to the lower wall and y = H to the upper wall, 
as shown in figure 1. Thus, the turbulent flow m the two-dimensional channel flow 
can be regarded as homogeneous in both streamwise (x) and spanwise (z) coordinates, 
while strong inhomogeneity is retained in the y-direction. For this case, one can 
define spectral functions for the first mode of velocity fluctuation, , and of 
pressure fluctuation, it , as follows. 

A 

<t, i (k 1 ,y,k 3 ,t) 


n(k 1 ,y,k 3 ,t) 

where k<| and k^ are wave numbers and i = /^T. The superscript indicating 
mode (1) is omitted in the spectral functions. Applying these definitions into 
eqs. (14) and (16), one obtains four complex equations with respect to the large 

A A 

eddy spectra, 4>- and ir. These equations can be further divided into eight equations 
for the real ancl imaginary parts defined according to the following notation. 


(2ir) £ 


(2tv) £ 


{/x<t>( 1 ^(x,y,z,t)}exp{-i(k^x + k^zHdx dz 

JT {/XTT^^x.yjZjtnexpI-U^x + k^zjldx dz 


(17) 
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(18) 


*2 = P 3 + iP 4 | 

*3 = P 5 + iP 6 

n = p 7 + iP 8 ^ 

The nonlinear terms require the convolution theorem (Lumley and Panofsky, 1964) 
during the transformation of the system of equations, eqs. ( 1 4 ) and (16), into the 
mixed, (k.| ,y , k^, t) , space. The spectral equations become 

A A A A 

ik-jir - (v + \> T )(<t>'j - k ct^) =0 

(19) 

- - (v + m t )(^ - k 2 J 2 ) = 0 

( 20 ) 

ikgU - (v + v T ) ( (}>” - k 2 <t>g) = 0 

( 21 ) 

ik ^ 4> ^ + $2 + lk 3^3 = 0 

( 22 ) 


« - - r a (D fu 

»-• " s "i S '! 


•ft • 


A A 


34 >r 

+ ik 1 U4, 2 + S 




34>q 

— + ik 1 U4» 3 + S 




where ( 


= d( )/dy 


and k 


k 2 + 


From equations (19) —(22 ) , we denote 



omd 

■i 


>} 


IT 2k* 1 'i i (k*’)J 1 (^ - k") 

* — CD 


d ♦,(£") « 

“d — *2 ( * “ k "> 


+ ik^4> 1 (S")i 3 (£c - jc")dk” dk^ 


(23) 


where k = (k^k^), k" = (k'j,k!p, and the double prime, ( )", when pertaining to the 
wave number variable denotes a dummy variable for the integration. 
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From a consideration of two-point velocity correlations and turbulent intensi- 
ties, a relation can be found between the spectra of the double velocity correla- 
tions and the 4>^ ' s as derived in appendix B. 

R ij )(k r y,k 3’ t ^^ + ^ = * i ( k 1 »y> k 3ft)4 , j(k' 1 \y, k g f t) (24) 

where 6 is the Dirac delta function and the superscript (1) indicates the contri- 
bution of the dominant mode. R^ is defined in 

R ij (r 1 » y * r 3» fc ) = // R i j( k 1 »y» k 3»t)exp{i(k 1 r 1 + k^Hdkj dk 3 (25) 

For a flow which is homogeneous in (x,z) -planes, 


R i j ( r 1 » y * r 3 » t ) = u i (x,y,z,t)u (x + r^y.z + r 3 ,t) 

= ]£ a^H^x.y.z.tH^U + r 1 ,y ,z + r 3 ,t) 
n=1 J J 

- a^^(x,y,z,t)4!^(x + r^y.z + r 3> t) (26) 

where r^ and r 3 are separation distances and = indicates simple truncation after 
the first mode. When r^ = ^ = 0 in eqs. ( 25) a nd (26), the two-point correlation 
reduces to the usual Reynolds stress tensor u^Uj(y,t). Thus, at time t, 

u i Uj (1) (y,t) = JJ 4> i (k 1 ,y,k 3 ,t)4>*(k 1 ,y,k 3 ,t)dk 1 dk 3 (27) 

where ( )* denotes the complex conjugate of ( ) . 

* 

The following structural quantities may also be calculated : (a) normal stress 

2 2 2 2 2 2 2 

intensities (u^/q ), where q = u + v + w , (b) shear stress intensity (uv/q ), 

(c) orientation of the principal axes of the large eddies (0), which is given by 
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0 = 1 tan 


( 28 ) 



2 2 

and (d) anisotropy (u /v ). 


CALCULATION PROCEDURES 


The objective here is to calculate the first mode of the velocity fluctuation 
and to investigate how much the first mode influences the turbulent intensities and 
the various other structural quantities described in the foregoing section. Atten- 
tion is paid to the skewness parameter, S, that appears in eq. (16), and its role in 
the solution for the first mode, 4>^ , and the structure of turbulence deduced 

therefrom. 

In order to solve eqs. ( 19)— ( 22 ) , after adopting a numerical algorithm one 
needs (a) the initial and boundary conditions, (b) the local mean velocity profile, 
and (c) the proper skewness factor as well as an eddy viscosity. In light of the 
difficulty in specifying the boundary condition for the pressure spectrum at the 
wall, the continuity relation is employed in the calculation over the Poisson's 
equation. Numerical results of the statistical quantities can then be compared with 
the measurements from Laufer (1951), for example, from which a particular flow 
condition is selected as 


Uo = 7.574 (m/sec) 
u* = 0.2891 (m/sec) 

H = 12.7 (cm) 

Re = UoH/v = 61600 

where Uo, u#, H, and Re are the mean velocity at the channel centerline, wall- 
friction velocity, channel width, and Reynolds number, respectively. 


Initial and Boundary Conditions 

A numerical solution for the large eddy spectra governed by the system of 
equations ( 19)— (22) is determined as an initial-boundary value problem m the (y,t)- 
space for various values of wave number, k. The initialization can be approximate 
because of the implicit nature of the algorithm and the goal of achieving a steady 
solution in the presence of a fixed mean strain. The initial conditions are chosen 
to possess reasonable spectral character and spatial distributions. The chosen 
initial distributions are as follows. 
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( 29 ) 


P 1 = Ay exp(-ky) > 
P 2 = Ay ( 0 . 5 - y)/k 



P E = A(1 - y)/k/B 
z> 

P 6 = P 5 ^ 

where A is a constant in the order of 1 , and the factor B in Pc has been 

2 2 2 3 

introduced so as to make the energy content of u , v , and w , that is, 

/ \ 

o 

isotropic initially. 

At the wall (y = 0), the no-slip condition requires 


P 


1 = 



(30) 


The boundary conditions for the pressure spectra (Py and Pg) are deduced from the 
v-component (normal to surface) equation applied at the wall (Moin, Reynolds, and 
Ferziger, 1978). The spatial derivative for the pressure fluctuation is then dis- 
cretized using a three-point one-sided formula beginning at y = 0. 


For the other set of boundary conditions, the flow field has been assumed to be 
symmetric with respect to the centerline (y = H/2), giving 


p . = p , = p , = = P . = p , = 0 (31) 

P 3 = P 4 = 0 (32) 

where ( )' denotes the derivative with respect to y. It should be pointed out that 
the flow field in the entire channel from y = 0 to y = H has been solved in a 
single instance with the no-slip condition imposed at both ends, y = 0 and y = H. 
The results show symmetric profiles for the u- and w-spectra, and antisymmetric 
profiles for the v-spectra (Pg and P^) with respect to the channel centerline. 

This justifies using the current boundary conditions at the centerline. 

The semi-implicit numerical scheme employed (Greenspan, 1974) utilizes a two- 
point backward differencing in time and a three-point central differencing in y. 

The nonlinear convolution integrals are treated explicitly by evaluating them at a 
previous time step when the solution is known. The numerical integration for these 
terms is carried out employing the trapezoidal rule over the wave number space, 
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(k'j,k"), at a given point, (k^,y,kg,t). This enables the formation of a system of 
matrix equations for all y at each advanced time level, where the coefficient 
matrix becomes block-tridiagonal and diagonally dominant. The inhomogeneous y 
coordinate is discretized as suggested by Murphy and Rubesm (1979) and the half of 
the channel is divided by 35 nonuniform grids. 

The wave number space, (k^kg), has also been divided into strongly nonuniform 
meshes. The wave number plane is covered with ( 17, 17) -grids where the values are 
equally spaced in the logarithmic scale for each wave number direction in the range 
between -10 and 10 (1/cm). It was found (from numerical experimentation) that wave 
numbers outside this range contribute so little energy to the first mode that they 
have negligible effect on the turbulent stresses. 


Mean Velocity Profile 


The mean velocity profile is approximated by a near wall Prandtl-Taylor model 
and a blending profile near the center plane. 

(a) U + = y + (y + < 12) 

(b) U + = 3.0 in y + + 5.5 (12 < y + < 760) 

(c) U/Uo = 1.0 + 0.068 log(y/d) (y + > 760) 

where U + = U/u*, y + = yu*/v, and d = H/2. The mean velocity profile in the outer 
layer, (c) above, has been introduced m this form for the purpose of matching with 
the law of the wall, (b) above, smoothly and, of course, with the experimental data 


Parameter Effects 

During the computation, a fairly small time step (normalized by the mean veloc- 
ity at the channel centerline and by the channel height) of about 0.001 has been 
used to ensure numerical stability and accuracy. Since the current system of spec- 
tral equations is nonlinear, an instability in any part of the solutions soon propa- 
gates into other solutions. Thus, if solutions for high wave numbers become 
unstable, even though solutions for low wave numbers are stable, the nonlinear 
integral causes the entire solution to grow infinitely as the iteration proceeds m 
time. 

Although 4> , governed by eqs. (19) -(22), generally can vary in time, the fully 
developed, steady mean flow and boundary conditions used here cause the solution to 
converge to a steady state which represents a stationary random field. We first 
studied the effect of various values of S on the solution when v^, = 0. Only a 
narrow range around the value of zero is found for the skewness parameter, namely 
| S | < 0.03, for which the solution does not grow rapidly. In the current computa- 
tion, the procedure has been continued up to 200 iterations in time to achieve an 
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accurate asymptotic solution. On a Cray operating system, the computation time took 
30 sec in the central processing unit for every iteration. 


NUMERICAL RESULTS 


First, this section will cover how the solution is affected by various values 
of the parameters, S and v T , in eqs. (19 ) — ( 22 ) . After a pair of those values is 
selected, numerical results obtained in the mixed, (k 1t y,k^,t), space are integrated 
over the wave number space to yield Reynolds stresses as a function of y. The 
computed stresses and structural quantities are compared with the experimental data 
of Laufer ( 1951 ) . 


Parameter Determination 

The skewness parameter, S, can be regarded as a parameter affecting the struc- 
ture of the solution, <j>. On the other hand, the primary role of the eddy viscosity, 
is to stabilize the growth of the solution, subject to production in a fixed 
mean velocity field. We pick a value of S first and then determine a correspond- 
ing value of which yields a steady-state solution. One may argue that a choice 

of a particular set of S and is not unique on the ground that other combina- 
tions of S and v T could also produce steady-state results. It has been found, 
however, that the solution is rather insensitive to the choices of the combinations 
of S and \> T that yield steady solutions (see below). 

Figure 2 shows the growth of the turbulent kinetic energy, integrated over the 
channel, for various values of S as a function of time in the absence of higher 
modes, = 0, and for a single case with = 18. With v = 0, the growth rate 
increases with increasing S, and it was found for |S| > 0.03 the solution grows 
so rapidly that it becomes unstable. For S = 0.01 and v„ = 18, however, the 
desired steady state m kinetic energy is achieved. Also, for S = 0.03 a steady 
solution for kinetic energy occurs when v^, = 22. 

In order to examine further whether each component of the kinetic energy has 
indeed reached a steady state for the above two sets of parameters, each component 
of the three turbulent intensities is integrated over the channel from y = 0 to 
y = d. The results are presented in figure 3 as a function of time. The results 
show that the u-component energy for S = 0.01 and v^, = 18 maintains a constant 
value, but for S = 0.03 and = 22 it continues to decrease slightly. For both 

cases the w-component continues to increase whereas the v-component decreases, 
again at a slow rate. The behavior of these different growth rates is believed to 
be attributed to the use of an isotropic eddy viscosity in eq. (15) and suggests the 
use of an anisotropic eddy viscosity or some other alternative (Hong and Murthy, 
1984a). Nevertheless, in view of the small growth rate m the v- and w-components, 
this behavior is believed to be relatively unimportant and no attempt was made to 
eliminate this continually varying anisotropy. The author favors the case of 
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S = 0.01 and \> T = 18 because it yielded a steady solution, more nearly steady 
than the other case. However, the calculated Reynolds stresses for these two sets 
of parameters were well within a few percent of each other as compared in figure 4. 
Thus, because the choice of which set to use is not critical to the results shown in 
the rest of the paper, S = 0.01 and v,p = 18 has been used. 


Eddy-Eddy Interactions 

The eddy-eddy interactions affecting the net production of the first mode will 
be shown as 
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The term "A" governs interaction with the mean flow. The term "B" corresponds to 
the transport of the '$> , whereas the term "C" represents the effects of the 

higher modes modeled by the J eddy viscosity. Figures 5 and 6 illustrate the nature 
of interactions among the eddies which are identified with "B" and "C" m eq. (33). 
For illustration, the large eddy-large eddy interaction in the P.j-equation (which 
is the real part of eq. (19)) is used, and is shown in figure 5 for three sets of 
wave numbers, 0.2, 0.5, and 1.0 (1/cm). It is found that the values of "B" in 
eq. (33) for wave numbers the same as, or less than, 0.2 (1/cm) are predominantly 
negative in the inner part of the boundary layer. A negative value of "B" refers to 
an energy supply; whereas a positive value indicates an energy dram because the 
value of Pi itself in the inner layer is negative for those wave numbers. The 
profiles of the nonlinear term show both types of behavior at the wave number about 
0.5 (1/cm) and positive behavior at higher wave numbers than 0.5 (1/cm). Thus the 
nonlinear eddy-eddy interactions for lower wave numbers cause energy gain, while 
those for higher wave numbers dissipate the energy. 

Shown in figure 6 are comparisons between "B" and "C" terms in eq. (33) for 
the Pi -equation. Figure 6(a) shows terms "B" and "C" at the low wave number 
ki =0.1 (1/cm) and figure 6(b) compares the same terms at the value of 
ki = 5.0 (1/cm), both for a fixed value of the wave number kg at 0.05 (1/cm). 
Profiles of the Pi are also provided m the two figures to indicate their behavior 
in y at the same wave numbers. An opposite sign in "B" or "C" from that of Pi 
implies energy loss, while the same sign implies an energy gain. The value of 
kg = 0.05 (1/cm) was chosen in these illustrations because it emphasizes the differ- 
ences occurring m the alternative ki's. Similar results are expected of "B" and 
"C" for other values of kg when ki is varied in the same manner. Figure 6(a) 
shows that at the same wave number (ki = 0.1) and y/d = 0.5 the turbulent energy 
transfer due to the first mode self-interactions is much smaller than that due to 
the rest of the modes when the latter are modeled with = constant. If other 
models were employed for the eddy-eddy interactions, this emphasis on energy loss 
might be reduced. The magnitude of the higher-mode interactions considerably 
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exceeds that of the first mode self-interactions in both figures 6(a) and 6(b) in 
order to dissipate the energy gained not only from the large eddy but also from the 
mean flow (term "A" in eq. (33)). Again from the signs of P^, "B," and "C" in 
figure 6, one can infer that the higher-mode interactions drain the energy consis- 
tently for all wave numbers, while large eddy self-interactions either supply energy 
for lower wave numbers or remove it for higher wave numbers. 


Turbulent Stresses 

The real part of the velocity spectrum for the u-component, P|, and that of 
the pressure spectrum, Py, is shown at two different locations of y, which is 
obtained at time t = 0.2. In figures 7 and 8, the P^ and Py-spectra are given in 
the first quadrant of the (k^ , kg) -space at fixed values of (a) y + = 7 and 
(b) y/d = 0.09. Generally, one can observe a smooth behavior in the distribution 
of Pi and Py. The value of Pi in figure 7 is greater at low wave numbers than at 
high wave numbers, and its distribution falls toward zero with the increasing wave 
numbers. As the value of y changes from y + = 7 to y/d = 0.09, the lower wave 
numbers contribute more energy to the Pi -spectrum in the inner layer region than in 
the viscous dominant region. In figure 8, the Py-spectrum, unlike P^ does not 
change its distribution much as y varies from one position to another. From these 
deterministic large eddy spectra, Pi through Pg, Reynolds stress components have 
been obtained as a function of (y/d) according to eq. (27). 

In figures 9-12, the Reynolds stresses obtained as a function of (y,t) are 
given at every 25 time steps to show the development in time. The solution adjusts 
itself quickly in time and the effect of the initial conditions appears to be 
minimal. In order to see whether the Reynolds stresses in figures 9-12 change their 
2 2 

profiles, the u - and v -profiles at t = 0.15, 0.175, and 0.2 are compared in 
figure 13 as a function of y/d. It shows that the u- and v-component intensities 
have indeed achieved equilibrium profiles for t > 0.15. The same observation has 
been made for other Reynolds stress components. 

For a detailed comparison, the Reynolds stresses are given as a function of 
y/d at time t = 0.2 in figure 14 along with experimental distributions taken from 
Laufer (1951). The first mode contributes approximately 30% of the observed inten- 
sities, although the shape agrees in general trend with the experimental distribu- 
tion. The use of isotropic viscosity has caused a spuriously higher proportion of 
the calculated w-component (than either the u- or v-component) in the contribution 
to energy. 


Structural Quantities 

The single-point structural quantities defined in Applications to Channel Flow 
were calculated from the Reynolds stresses and are compared with corresponding 
quantities obtained from the measurements of Laufer (1951). In spite of low 
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intensity levels of the large eddy as shown in figure 14, the normalized structural 
quantities m figure 15 agree well with the corresponding experimental quantities 
except for the anisotropy of the flow in the outer part of the boundary layer shown 
m figure 15(d). It is believed that this latter result is the consequence of the 
imposed symmetry at the channel centerline, that is, v = 0 . 


LARGE EDDY STRUCTURE 


An important quantity in the study of turbulence structure is the two-pomt 
velocity correlation which depends on the distance separating the velocities at two 
different positions. Other useful quantities include the fluctuating velocity 
vector and stream functions from which a flow pattern may be constructed. We shall 
present these quantities as deduced from the large eddy computation and discuss 
their implications in this section. 


Two-Pomt Velocity Correlation 

To demonstrate the character of the velocity correlations that can be evaluated 
from the computation of the first mode, four correlations, R^, R 22 , R 33 , and R 12 , 
are presented here. The points are separated by r-|, r 2 , and r^ m the x-, y-, 
and z-directions, respectively, at a fixed time, t = 0.2. First, we have calculated 
correlations between two velocity components which are separated by (r 1 and r^) m 
the homogeneous plane (y = constant). Then consideration is given to correlations 
between velocities separated by r 2 m the y direction. 

In figures 16 through 23, the two-point velocity correlations (obtained from 
eq. (25) and normalized by the maximum value at each y location) are shown at two 
planes across the channel, namely y + = 7 and y/d = 0.5, as a function of separa- 
tion distance (r^r^). At y + = 7 m the viscous sublayer the correlations fall 
off rapidly to a zero value, but oscillate thereafter as the separation distance 
increases in both x- and z-directions . A negative region in the correlation sur- 
face indicates that the velocity fluctuation in that region is generally opposite 
from the one at the origin. Thus, the distance between the two adjacent peaks m 
figures 16-19 reflects a macro-length scale of motion in the region. On the other 
hand, the correlations at y/d = 0.5 drops slowly over the distance of approxi- 
mately a channel half-width. In general, one can infer that the normal velocities, 
v, have shorter correlation lengths than the other velocity components. When the 
correlations in figures 16-23 are compared with the ones observed experimentally 
(Comte-Bellot, 1963), the distance for which the correlation changes its sign or 
reaches its minimum is longer in the current calculations than m the measurements. 
In general, these correlations exhibit the different character among the components 
in R x j and between different locations m y. For example, the longitudinal 
correlation, R^ in figure 20 , and the lateral correlation, R 22 m figure 21 , 
display the typical behavior of these quantities; that is, the R^ asymptotes to a 
zero value while the R 22 crosses the zero line in a short distance and comes back 
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to R 22 = 0. But all correlations are symmetric with respect to r 1 and ro for 
all y's, which is consistent with the homogeneity assumption of the flow in the 
(x,z) -plane. 

An interesting way of examining the correlation coefficients is to introduce a 
length scale, 2-, deducible from the correlations in the form 

*i (y) = / / R ii (r 1’ y,r 3 )dr 1 dr 3 (34) 

*0 # o JO 

Here the subscript, i, may be 1, 2, or 3 corresponding to the x-, y-, or 
z-coordinates , respectively. The length scales obtained from eq. (34) for different 
" 1 " are compared across the channel in figure 24. One observes the distinct anisot- 
ropy in the length scales among the directional components. Also the magnitude of 
the length scales is on the order of the channel half-width, d, with the largest 
x-component and the least y-component scales. The length scales defined this way 
represent the size of the large eddy (Townsend, 1976, p. 50) along y, where the 
size in the ith direction is proportional to each component of the length scale, 

2 . 

1 

Next, normal correlations of the velocities are examined as a function of 
(y,r 2 ) m the form 


u. (y)u (y + r p ) 

= = ( 35 ) 

{u ?«a* 

In the calculation, the R li are evaluated at t = 0.2 for every (x,z)-point in a 
range which covers about two periods m each of the x- and z-directions. Then the 
correlations are averaged over the (x,z)-plane. Such correlations are shown at two 
values of y in figures 25 and 26. At y + = 7 in figure 25, the u- and 
v-correlations approach asymptotically a near zero value within a distance of about 
O.ld, whereas the w-correlation drops to a zero value within 0.05d. On the other 
hand, the correlations at y/d = 0.5 decrease mildly with increasing separation 
distance as shown in figure 26. The shapes of the three correlations m figure 26 
are also very similar to one another except in the region very close to the wall. 

The anisotropy among the three correlations is more severe at y + = 7 than at 
y/d = 0.5, but the correlation is weaker m the viscous region than m the outer 
layer. It means that the velocities are more likely to be in the opposite sign in 
the viscous region than in the outer layer. 

Again length scales are assigned to the correlation coefficients defined in 
eq. (35) m the following form. 



t, ll (y,r 2 >dr 2 


(36) 
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The length scales, L^, then represent the degree of mixing of the fluid in the y 
direction associated with each velocity component. Figure 27 shows the length 
scales obtained from eq. (36) as a function of y/d. On the same figure, a mixing 
length scale deduced from 

™ * -»*« (§f <37 

is provided for a qualitative comparison between L x and £ m and is denoted by 
points (•). The values of the scale, J, , are taken from Goldstein (1950, p. 357) 
for a pipe flow at the same Reynolds number as in the current channel flow. Thus 
the length scales, L i? are comparable to the conventional mixing length scale in 
magnitude as well as in distribution. But, as can be observed in figure 27, the 
x-component length scale, L i? becomes locally negative m the viscous region. This 
is perhaps because the effect of small-scale turbulence is excluded in that region 
where the influence of small eddies is comparable to that of large eddies. 

As an alternative to the integral scale, L^, as defined in eq. (36), one can 
define a length scale as the distance for which the correlation crosses a zero 
line. The slope of such a length scale may be found to be steeper (Glushko, 1965) 
than the current slope of about 0.4 in the inner region (y/d < 0.2) m figure 27. 


Typical Velocity Fluctuations 


Characteristics of the turbulent velocity field can be displayed, in the 
ensemble-averaged sense, from the first mode, 
field corresponding only to the first mode as 
from 


If one writes at typical velocity 
' then u. ' may be obtained 


u. 

i 


u[ 1 }(x,y,z,t) = y X ^ ^4>^(x,y,z,t) 


= JJ 4> l (k 1 ,y,k 3 ,t)exp{i(k 1 x + k^zHd^ dk^ (38) 

Contour plots of the typical velocity fluctuations can then be drawn m a physical 
domain between -0.5H and 0.5H for both x- and z-coordinates . To show the turbulent 
motion in the near wall region, the contour plots of the three velocity components, 
U (D, v' 1 ', and w^, are presented at y + = 7 in figures 28-30. In the figures, 
solid lines indicate positive values and dashed lines indicate negative velocities. 
The inner curves possess higher values than outer curves enclosing them, and the 
increment between contour curves is uniform throughout the figures. The solid and 
dashed lines appear alternatively in both x- and z-directions showing a periodic 
occurrence of the "thumbprints." The contour curves are more circular than similar 
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results from the direct simulation of the flow (Moin and Kim, 1982) where the con- 
tour lines form elongated curves in the streamwise direction. The contour plots 
exhibit periodic formations in both directions and the spatial periodicity is shown 
to be one-half of the channel height (H/2), although there are also weak secondary 
periods within the principal period. This is a consequence of the Fourier transform 
utilized in the calculation procedure. One advantage of the spectral method over 
the finite difference method is that the latter requires a period length in the 
homogeneous directions when imposing the boundary conditions while the former does 
not need it. The magnitude of the period in the spectral method is determined 
implicitly as a function of the large eddy spectral content and is dependent on the 
nature of turbulence rather than artificial geometry. Since the flow domain m the 
figures covers two periodic lengths in both directions, it is considered adequate in 
studying any feature deduced from the region as a typical property of the flow as a 
whole. 


Stream Function 

The velocity vector of the first mode is plotted in figure 31 m the 
(y,z)-plane at two values of x, that is, x/H = 0.3 and 0.4. The choice of 0.3 and 
0.4 is not significant since the turbulent motion is shown to be periodic and repet- 
itive m the previous section. The two locations are chosen simply because they 
seem to possess some character of the flow. The flow at 0.3 and 0.4 appears to be 
chaotic m figure 31. In order to grasp the flow field associated with the velocity 
vectors, a stream function, 4^, is defined m the (y,z)-plane as follows. 


^ 1 ( x , y , z ) 




(x,y,z)dy 


(39) 


The stream function, ^ , represents the turbulent flux of mass between the surface 
and a point in space. A positive value implies the bulk of the flow moving in the 
negative z-direction and a negative value implies movement in the positive 
y-direction. The contours of the stream function are drawn m figure 32 and, in a 
three-dimensional format, in figure 33, where the geometrical extent of the flow 
field in the (y,z)-plant is the same as m figure 31. The solid lines bear positive 
values, and the dashed lines correspond to negative stream functions. Figure 32 
shows contour curves of 4^ obtained at (a) x/H = 0.3 and (b) 0.34. The stream 
function contours when chosen at 0.3 and 0.34, instead of 0.3 and 0.4, reveal 
clearly the turbulent motion in the spanwise direction as well as that in the 
streamwise direction. First, the bulk motions of turbulent flow push each other m 
the spanwise direction, resulting in a counter-rotating flow. The centers of each 
bulk motion are located at the regions of highest contour values in figure 32. 
Secondly, the line connecting points A-E in figure 32 rises to the line A'-E' in the 
adjacent downstream plane. The values of the stream function remain unchanged as 
the flow moves from A to A' , and B to B', etc. The streamlines connecting the two 
lines are elevated by approximately 0.06 m y/H between 0.3 and 0.34 in x/H. 
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Thus one can estimate that the large eddy structure is inclined with respect to 
free-streamline by an angle of arctan( 0.06/0.04) = 56 (degrees). 

In figure 33, the levels of the stream function contours remain unchanged as 
the x/H progresses from 0.3 to 0.5. The streamlines at x/H = 0.5 have been 
added primarily to show the periodic nature of the motion in the homogeneous direc- 
tion. The main feature in figure 33 is the lifting of the streamlines as the flow 
moves on from x/H = 0.3 to 0.4. If one compares figure 31 with 33, the flow in the 
inner region (y/d < 0.3) at x/H = 0.3 is more energetic than at x/H = 0.4. This 
is because the flow has gushed up toward the outer region by a movement of counter- 
rotating vortices at x/H = 0.3 and as a result has lost its energy by the time it 
reaches x/H = 0.4. The lower portion of the flow field m figure 31 is indeed very 
quiet and represents a very weak motion. After resting briefly at around 
x/H = 0.4, the flow is then ready again for another leap at around x/H = 0.5. 

When attention is focused on the streamlines in the viscous region, as in 
figure 34 where physically the spanwise extent is six times the vertical extent, two 
opposing streamline patterns emerge. The distance between the centers of the domi- 
nant motions is about 140 wall units (normalized by the wall-friction velocity and 
the molecular viscosity). This is higher than the 50 wall units estimated by 
Blackwelder and Eckelmann (1979) for the spanwise distance between a pair of 
counter-rotating streamwise vortices. Nevertheless, the 4>. shown on figure 34 
again demonstrates the existence of a counter-rotating vortex motion which intrudes 
the viscous layer. These are the same motions which have been advocated by Bakewell 
and Lumley (1967) and by Moin (1984). 

In addition to the previously defined stream function, i^, in the (y,z)-plane, 
a second stream function, tl^, may be considered in the (x,y)-plane to attempt to 
explain the manner in which the large eddy behaves in the streamwise direction. 

Now, if the second stream function, is expressed in the form 


♦ 2 < x *y» z ) 




(x,y,z)dy 


(40) 


then the ^ describes the flux of the turbulent flow in the streamwise direction 
across a vertical line from the wall (y = 0) to a point of interest (y). The curves 
pushed toward the positive y-direction imply that the flux of turbulent flow is 
going into the positive x-direction, while the curves dipped toward the bottom wall 
indicate that the turbulent flow is returning back against the mean flow even though 
the total velocity is still positive. In figure 35, if one examines the flow region 
between x/H = 0.3 and x/H = 0.4 only, the streamlines bounded by the first two 
(x,y) -planes from the low z values are indented toward the bottom. On the con- 
trary, the streamlines surrounded by the last two (x,y)-planes show the opposite. 

One can infer, therefore, that the motion of the large eddies pushes the fluid 
forward, stagnates, and then pulls it backward, forming a shape analogous to a 
horseshoe or hairpin (Willmarth and Tu, 1967). However, the structure inferred from 
figures 33 and 35 covers a greater spatial extent than the one suggested, for 
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example, by Hinze (1975), and virtually the whole flow field becomes affected where 
such structure arises. The structure occurs sporadically in space, with a principal 
period equaling the channel half-width in the homogeneous directions. Within the 
principal period smaller, sub-harmonic motions are also observed. 


LARGE EDDY STATISTICS 


While correlation functions, vector plots, or stream functions of the first 
mode reveal an essential nature of the turbulenty motion, another side of turbulence 
can be exposed by vorticity fluctuations and probability functions. To enhance our 
understanding of the large eddies, we further investigate the velocity fluctuations 
from different points of view and their respective statistical properties. We shall 
briefly examine the hodograph transformation of the large eddy velocity components. 


Large Eddy Velocity Field 

The same velocity fluctuations, obtained from eq. (38) and discussed m Typical 
Velocity Fluctuations, are then displayed in the (y,z)-plane mainly to show their 
variations in the inhomogeneous space, y, for successive locations of z. Most 
fluctuations are concentrated m the wall region as shown in figures 36-38 where the 
energy production is large. One notable feature is how the fluctuation at the wall 
is related to the outer edge of the boundary layer, which is reminiscent of the 
ejection or inrush (Hinze, 1975, p. 666). The ridges m figure 37 may be identified 
with the ejection originating at the top of each hill, and the valleys can be 
related to the inrush. The two phenomena occur intermittently in the horizontal 
(x,z)-plane but are connected throughout the y-direction. This shows that the 
turbulent transport transmits most of the energy or the momentum through the trans- 
mittal routes m the normal direction. 


Vorticity Fluctuation 

Although an eddy is a different concept from a vortex, efforts have been made 
to trace the eddy from the dynamics of the vortex motion. The vorticity is further 
distinguished from the vortex but can be indicative of turbulent motion, especially 
through the formation of a vortex tube. The streamwise component of the vorticity 
fluctuation is calculated from the large eddy and is shown in figure 39. For quali- 
tative study, however, the differentiation of the velocity is performed in the 
physical space using a central differencing. If one wishes higher accuracy, one can 
differentiate the velocity spectrum analytically and then integrate it over the wave 
number space. The streamwise vorticity reveals highly concentrated counter-rotating 
vortices along the spanwise direction. The vorticity in the outer layer is surpris- 
ingly quiet m contrast to the energetic velocity fluctuations m that region. 
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Hodograph Transformation 


With the anticipation that the hodograph transformation may hold a clue for the 
description of the large eddies, the velocity flow field of the first mode in a 
finite domain in the physical space has been transformed into the so-called hodo- 
graph-(u, v,w)-space. A flow volume taken from the (x,y ,z)-space at a fixed time, 
t = 0.2, as shown in figure MO, is represented by 35 grid points along each of 
36 lines sampled. Velocity components at each grid point along the vertical lines 
are then mapped into the hodograph space, as shown in figure Ml, and the points are 
connected m the direction of the arrows. The result is 36 tangled lines m the 
(u,v,w)-space, forming a finite region in the new space. Each trajectory in the 
hodograph space consists mostly of loops with relatively large radii of curvature. 

As noted by Deissler (198M) m the description of a flow with all sizes of eddies 
transformed into a phase space, the trajectory for the small-scale turbulence might 
have included loops with smaller radii of curvature than the ones shown m fig- 
ure Ml. It is also noticeable that there are no cusps m the trajectories due to 
inherent large scales associated with the first mode. In the figure, all the paths 
are heavily concentrated around the zero-velocity point, and the trajectories follow 
the curve, u = -v, when projected on the (u,v)-plane as shown in figure M2. 

By the same way, streamlines passing through the (x,y)-plane at z = 0 are 
projected into the hodograph space. The grid points from which velocities are taken 
are shown in figure M3, and the velocities are then mapped into the hodograph space 
m figure MM. The transformed shape in figure MM is obtained after connecting all 
the points in both x- and y-directions . The trajectories again show very much the 
same shape as found earlier in figure Ml from another volume of flow. It is evident 
from the two figures that any streamlines in the turbulent channel flow, when con- 
verted into the hodograph space, would lie on the surface enclosed by the 36 arrowed 
lines in figure Ml. Thus the large eddy velocity components in the channel flow 
form an ellipsoid in the (u, v,w)-space whose major axis is aligned with a diagonal 
represented by u = -v = w. This is consistent with the proposed description of a 
large eddy by Townsend (1956) and other experimental observations in which correla- 
tions between u and v, for example, are most likely negative leading to negative 
shear stress. 


Probability Distribution Function 

In previous hodograph figures as well as velocity contour plots, it is apparent 
that the velocity fluctuations pass through the zero point (u = v = w = 0) more 
often than any other points. When 2601 (= 51 x 51) grid points were sampled in an 
(x,z)-plane at y + = 7 or y/d = 0.5, the number of points at which the velocity 
falls into a certain threshold between maximum and minimum velocities shows a typi- 
cal near-Gaussian distribution as in figure M5. The distribution in the homogeneous 
plane slightly deviates from the Gaussian and the degree of departure from it indi- 
cates the degree of skewness of turbulence. The skewness is believed to be a result 
of the nonlinear transport term m eq. (16) which was necessary for the energy 
balance of the large eddies. In figure M5, two propability distribution function 
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(PDF) curves are presented for two different locations of y, and the skewness 
changes its sign as y increases. This agrees with experimental trends (Kreplin 
and Eckelmann, 1979). When the u-component velocity is sampled along a vertical 
line in y at two locations of (x,z) = (0,0) and (0. 1,0.1) as illustrated in fig- 
ure 46, the resulting PDF is not even close to a Gaussian distribution, thereby 
confirming that the PDF in the inhomogeneous space is not normal as expected. 

The statistical results discussed m this section are obtained from the first 
mode; nonetheless they disclose salient features of the entire turbulent motion for 
a given flow, either experimentally observed or derived from analysis. It is thus 
certain that the first mode retains most of structural and statistical information 
embedded in an actual flow. 


TURBULENCE CORRELATIONS 


The nature of the orthogonal decomposition of the velocity fluctuation into 
random and deterministic variables in the Large Eddy Interaction Model could enable 
one to calculate various turbulent correlations directly from the solution of ' 
and ir' n . As an illustration, we now consider two important processes appearing in 
the Reynolds stress equations; that is, the pressure-strain rate correlation, tt , 
and the rate of dissipation of the turbulent kinetic energy, e^, where 1 " 



We recall the expansion in the form 
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. . . , are the solutions of the coupled dynamical equations 
(see eqs. (10) and (13)), but are not necessarily mutually 
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in the form 


Due to the randomness in the a 's, one can express it., and e., 
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it i j = Lc^ir j + a 2 it S^j + ...J 


where 




e ij 


n=1 


2v E 

n=1 




3X, 


3X, 


(45) 


(46) 


.(n) 


3 *! n) 

l 

3X . 


3 «J> 


(n) 


3Xi 


Concurrently with the homogeneity assumption in the channel flow case, ir. , and e . 

* ^ IJ IJ 

can also be expressed as a function of spectra, 4> . and it, as well (Hong and Murthy, 
1984b). 1 


Therefore the Large Eddy Interaction Model can be utilized as a guide to test 
and to aid in closing various turbulence models. 


CONCLUDING REMARKS 


The current paper deals only with the first mode and the balance of a large 
eddy regarding its origin, maintenance, and destruction . Therefore, the question 
still remains untouched as to the number of modes in the series required to repre- 
sent the turbulence field sufficiently well to predict engineering quantities. The 
first mode is shown, however, to be so significant that it supplies about 30$ turbu- 
lent kinetic energy and possesses a structural and statistical character that 
matches well with the experimental trends of the entire turbulence field. 

When the turbulent transport is truncated to only the interactions between the 
components of the first mode (v^ =0), it is found that the energy drain through the 
first mode is insufficient to balance the energy gam from the mean motion. The 
kinetic energy of the first mode then grows monotonically without bound. The non- 
linear interactions of higher modes are thus necessary in the dynamics of the first 
mode. However, the isotropic eddy viscosity used in the current closure yields an 
anisotropic growth rate, though small, among the three normal intensities. To 
improve the closure with respect to the higher mode-interactions, an alternative 
approach may be to solve the first and second modes simultaneously and to model the 
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interactions resulting from the rest of the modes. In that case, the different 
roles of the first two modes in forming a complete flow should become clear and help 
answer, in part, the earlier question of the number of modes essential for the flow. 

It is briefly mentioned in the last section that the Large Eddy Interaction 
Model can provide an independent framework to check models for various turbulent 
processes. The soundness of such a test lies in the fact that the Large Eddy Inter- 
action Model requires a closure assumption for only the nonlinear eddy-eddy interac- 
tion terms, and the empirical parameters are chosen on the basis of physical consid- 
eration of the large eddies. Now that the first mode is shown to represent the 
structural character of turbulence, it remains to be determined whether the first 
two modes are adequate to capture most of the turbulence energy. 
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APPENDIX A 


ALTERNATIVE EXPANSIONS FOR VELOCITY FLUCTUATION 


If the instantaneous velocity field is decomposed into a time-mean value and 
its fluctuation, then 


0 . 

1 

= U + u. 

l l 

(A.1) 

<0.> 

1 

= U. 

l 

(A. 2) 

<u.> 

= 0 

(A. 3) 


where < > denotes a time average. The velocity fluctuation can be expanded in 
terms of orthonormal functions as opposed to eq. (2). 

00 

u i (x,t) = £ c* n (t)<t/ n *(x) (A. 4) 


Notice that the orthonormal functions depend on space. The random coefficients are 
functions of time and are subject to <a (t)> = 0 by definition. Then the expan- 
sion, eq. (A. 4), satisfies the premise or eq. (A. 3), that is 

<u . > = y. o (t)^ n) (j) 

i in i 

n= 1 

= E < a n (t)>4>^ n) (x) 

n=1 

= 0 (A. 5) 

Thus, the expansion m eq. (A. 4) is consistent with the definition of the velocity 
fluctuation m eq. (A. 3). 

Similarly, if U^ in eq. (A.1) is a spatially averaged value, then 

{U i > = U 1 (A. 6) 

{u x } = 0 (A. 7) 

where { } represents a space average. The velocity fluctuation then can be 
expanded in the form 
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with 


{a n (£)> = 0 (A. 8) 

so that the expansion in eq. (A. 8) satisfies the definition of u^ 

( u 1 > = |E a n (x)<J>[ n ^(t)| 

= E {a n (x)}4»j n) (t) 

n=1 n 1 

= 0 (A. 9) 

Among the three types of expansions for u^, which are eqs. (2), (A.4), and 
(A. 8), eq. (2) is the most "fundamental" way of expanding the velocity fluctuations. 



APPENDIX B 


RELATION BETWEEN VELOCITY SPECTRUM AND SPECTRUM OF TWO-POINT VELOCITY CORRELATION 


If mean flow is homogeneous in both streamwise (x) and spanwise (z) coordi- 
nates, as is the case of two-dimensional channel flow, one can define spectra of the 
first mode, 4>^ , and of the two-point spatial velocity correlation as follows. 

= jhtSL A4>[ 1) (x,y,z)exp{-i(k 1 x + k^z)}dx dz (B.1) 


R ij (k 1 ,y , k 3 ) = - ^ JJ R i j(r 1 ,y,r 3 )exp{-i(k 1 x + k 3 z)}dx dz (B.2) 


or, inversely, 


/X4>[^(x,y,z) = JJ <J> i (k 1 .y.k^exptM^x + k 3 z)}dk 1 dk 3 (B.3) 


R ij( r 1 » y » r 3 ) = JJ R iJ (k 1 ,y,k 3 )exp{i(k 1 x + k 3 z)}dk 1 dk 3 (B.4) 


Here all variables associated with 4 k and R^j depend on time implicitly. 
Now, consider a product of 4 k and 4 >j in the form 


4 > i (k 1 »y , k 3 )4»j(k”,y,k”) 


^ - 1 JJ 1 ^(x,y,z)exp[-i(k 1 x + k 3 z)]dx dzj- 

x i — ff /x4k ^(x",y,z")exp[-i(k' 1 'x" + k"z")]dx" 



JlfXT^ X <’[ 1 )(x >y> z) *j 1 )(x ">y> z,,) 


x exp[-i(k,|X + k 3 z + k'jx" + k^z")]dx dz dx" dz" 


(B.5) 
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Let's set x" and z" as 


x" = x + r. 


z" = z + r. 


(B.6) 

(B.7) 


then 


* a. ( k i , y » k 3 ) * j ( k i , y » k 3 ) 


(2tt)‘ 


f 


.} 


(2tt) 


2 R x j ( k” ,y , k” ) JJ exp{-i[(k 1 + k'j)x + (k^ + k^)z]}dx dz 


(B.9) 


x exp[-i(k'jr 1 + k^r^)]dr 1 dr^ * exp[-i(k.j + k" )x - i(k^ + k^)z]dx dz (B.8) 

The integrand in the last expression can be replaced by the two-point velocity 
correlation, Rjj, by way of 

R ij(r t , y , r 3 ) = X^^^U.y.zH^U + r r y,z + r 3 > 

yielding 

4> 1 (k 1 ,y,k 3 )4. J (k , 1 ',y,k^) 

= X[.{^7 JL V r i' y ' r 3> ' exp[ - :<k i r i * lt 3 r 3 )ld '-i dr 3 

A 

x exp[-i[(k 1 + k!J)x + (k 3 + k^zjjdx dz 
1 


(B. 10) 


Fourier transform of a constant is then the Dirac delta function, that is 

6(& + &") = — ff exp{-i[(k. + k")x + (k_ + k")z]}dx dz (B.11) 

(2 *r id 


or 


1 



6(k + &")exp{{i< + k") 


x}d(k + &") 
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By substituting eq. (B.11) into eq. (B.10), one arrives at 


^> i ( k 1 »y» k 3)'l>j(k' 1 ',y,k^) = R iJ (k' 1 ',y,k^)s(^ + £") 

or 


R iJ (-k 1 ,y,-k 3 ) = <t> i (k 1 ,y,k 3 )<o j (-k 1 ,y,-k 3 ) (B.12) 

which is given in eq. (24) and leads subsequently to the following via eqs. (25) 
and (26). 


u i u j( y ) = JJ 4> x ( k ^ , y , k 3 ) j C ~ k i ,y,-k 3 )dk 1 dk- 


(B. 13) 


Here one can utilize the following symmetry condition implied in eqs. (B.3) and 
(B.4) for the integration. 


*y»~ k 3 ) 

A A 

4> 1 (-k r y, k 3 ) = <t>*(k r y,-k 3 ) 

R iJ ( k 1 »y» k 3 ) = R* J (- k 1 ,y,~ k 3 ) 

That is, the real parts of <J> and R are even functions, and the imaginary parts 
are odd functions with respect to the^wave number, 

Another way of obtaining u u (y) may be to (a) calculate <t)^(x,y,z) and 

( 1 ) 1 J 1 

4>, (x,y,z) separately from eq. (B.1), (b) multiply them directly in the (x,z)-plane 

J ( i ) ( 1 ) 

for <t>: 4> , and (c) take a spatial average over a period which the product would 

possess in J both x- and z-coordinates. Thus, the alternative will be more tedious 
than the analytical methods in the above, eq. (B.13), and the solution depends on 
the averaging technique which is unnecessary in the spectral method. 
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Figure 1.- Coordinate system for a two-dimensional channel flow. 


Figure 2 



- Turbulent kinetic energy vs t (= tUQ/H) for various values of skewness 

parameter, S. 





Figure 3.- Turbulent energy growth in time among three components for two sets of 
parameters: solid line for (S = 0.01, v, = 18); dashed line for (S = 0.03, 

v T = 22).. 
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S= 0.01, Vj= 18 
S = 0 03, v T = 22 



Figure 4.- Comparison of Reynolds stresses calculated, based on two sets of param- 
eters: S = 0.01 and = 18 (solid line); S = 0.03 and v T = 22 (000). 


O kl =k3 = 02 1/cm 
V kl =k3 = 05 1/cm 
□ kl = k3 = 1 0 1/cm 


2 4 6 

y/d 


8 1 0 


Figure 5.- Profiles of S(3/3x.) ^)(0<^)> versus y/d for various wave numbers 
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Figure 6(a).- Profiles of large eddy self-interactions and modeled higher mode 
interactions versus y/d for ki =0.1 and ko = 0.05 (1/cm) at t = 0.2. 
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Figure 



2 2 

11.- Development of spanwise turbulent intensity {w /U ) in the (y,t)-space. 



2 

Figure 12.- Development of shear stress (uv/U Q ) in the (y,t)-spaee. 
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I V 2 profiles at t = 0.15 (o), 0.175 (*), 

0.2 (a). 
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0, deg -uv/q' 


Figure 15.- Structural quantities calculated from the first mode versus y/d at 
t = 0.2: ••Laufer (1951); (a) uv/q 2 , (b) uv/l/u 2 ^ 2 , (c) 0 in equation (28), 


and (d) v^/u 2 . 
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Figure 16.- Two-point velocity correlation, obtained from eq. (25), as a 

function of (ri,rq) at y + = 7. 






Figure 17.- Two-point velocity correlation, obtained from eq. (25), as a 

function of (r^r^) at y + = 7, 





Figure 19-- Two-pomt velocity correlation, R 12 obtained from eq. (25), as a 

function of (r^,r^) at y + = 7. 
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Figure 20.- Two-pomt velocity correlation, obtained from eq. (25), as a 

function of (r^r^) at y/d = 0.5. 


Figure 21.- Two-point velocity correlation, obtained from eq. (25), as a 

function of (r^r^) at y/d = 0.5. 
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Figure 22.- Two-pomt velocity correlation, R 33 obtained from eq. (25), as a 

function of (r^,ro) at y/d = 0.5. 



Figure 23.- Two-point velocity correlation, obtained from eq. (25), as a 

function of (r^,r^) at y/d = 0.5. 



Figure 24.- Integral length scale of large eddy defined by eq. (34) vs. y/d. 
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x x-COMPONENT (R 1t ) 
V y-COMPONENT <R 22 ) 
O z-COMPONENT (R 33 ) 



Figure 25.- Two-point velocity correlations obtained from eq. (35) as a function 

of r 2 at y + = 7. 



Figure 26.- Two-pomt velocity correlations obtained from eq. (35) as a function 

of r 2 at y/d = 0.5. 
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Figure 27.- Integral length scale of large eddy defined by eq. (36) vs. 

••Goldstein (1950). 



Figure 28.- Contour plot of u^ in the (x,z)-plane at 


y + = 7. 


y/d: 
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Figure 30.- Contour plot of v-/ 1 ^ in the (x.z)-plane at y 




Figure 31 



Velocity vector of large eddy in the (y,z)-plane at (a) x/H = 0.3 and 

(b) x/H = 0.4. 
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Figure 33.- Stream function ^ (eq. (39)) contours in the (y f z)-plane at 

x/H = 0.3, 0.4, and 0.5. 
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Figure 34.- Stream function il>. contours m the (y,z)-plane in the viscous region 





Figure 35.- Stream function 


\|>2 (eq. (40)) contours in the (x,y) -plane at four 
locations of z. 



Figure 36.- Large eddy velocity fluctuation, u^, in the (y,z)-plane at 

x/H = 0.12. 
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in the (y,z)-plane at 


Figure 37.- Large eddy velocity fluctuation, v v ', 

x/H = 0.12. 



Figure 38.- Large eddy velocity fluctuation, w^, in the (y,z)-plane at 

x/H = 
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Figure 41.- Traces of velocities along the lines shown in figure 40 at t = 0.2. 
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Figure 42.- Projection of velocity traces in figure 41 onto the (u,v)-plane, at 

t = 0.2. 
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Figure 43.- A plane in the physical space at z = 0. 
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150 

Figure 44.- Traces of velocities taken from every grid point in figure 43, at 

t = 0.2. 
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Figure 45 




(b) U (D 


Probability distribution function of the u-component of a large eddy 
in the homogeneous space: (a) y + = 7, (b) y/d = 0.5. 
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Figure 46.- Probability distribution function of the u-component of a large eddy 
along y at (a) x = z = 0 and (b) x = z = 0.1H. 
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